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ABSTRACT 

In order to investigate whether massive stars form similarly to their low-mass coun- 
terparts, we have used the standard envelope plus disc geometry successfully applied 
to low-mass protostars to model the near-IR to sub-millimetre SED and several mid- 
IR images of the embedded massive star IRAS 20126+4104. We have used a Monte 
Carlo radiative transfer dust code to model the continuum absorption, emission and 
scattering through two azimuthally symmetric dust geometries, the first consisting 
of a rotationally flattened envelope with outflow cavities, and the second which also 
includes a flared accretion disc. Our results show that the envelope plus disc model 
reproduces the observed SED and images more accurately than the model without a 
disc, although the latter model more closely reproduces the morphology of the mid-IR 
emission within a radius of 1.1" or ~1800au. We have put forward several possible 
causes of this discontinuity, including inner truncation of the disc due to stellar irra- 
diation, or precession of the outflow cavity. Our best fitting envelope plus disc model 
has a disc radius of 9200 au. We find that it is unlikely that the outer regions of such 
a disc would be in hydrostatic or centrifugal equilibrium, however we calculate that 
the temperatures within the disc would keep it stable to fragmentation. 

Key words: accretion, accretion discs - radiative transfer - stars: formation - in- 
frared: stars - ISM: jets and outflows - reflection nebulae 



1 INTRODUCTION 

Does the formation of a massive star differ significantly from 
that of a low mass star? This is the question at the core of 
all studies of high mass star formation to date. The fact 
that massive stars can have Kelvin-Helmholtz time-scales 
shorter than their accretion time-scales ( |Shu et al.| [l987), 
leads to the consequence that their main-sequence radiation 
can affect their own formation. Hence, several processes such 
as radiation pressure (e.g. Yorke||2u02 ) and ionization (e.g. 
Keto 2002) may halt, decrease or alter accretion on to the 
star. 

With the hypothesis that these differences do not radi- 
cally change the formation of massive stars, our main aim is 
to find observational evidence to suggest the contrary. In this 
paper, we aim to determine whether the standard envelope 
plus disc model used to reproduce the spectral energy dis- 
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tributions (SEDs) of low mass stars (e.g. Adams et al.|198^ 



Kenyon et aT|l993| |Whitney et al.|2003| |Wood et al.|2 002 ) 

is able to reproduce the SED and observed infrared images 
of the nearby forming early B-type star IRAS 20126+4104. 



Keto fe Zhang| fl2010[ henceforth KZ10) have re- 
cently modelled the line profiles of several millimetre- 
wavelength transitions of NH3, C 34 S and CH3CN toward 
IRAS 20126+4104, fitting these data with two models, the 
first consisting of an azimuthally symmetric rotationally flat- 
tened envelope, and the second also including a flared Ke- 
plerian disc. They found that the addition of a disc was 
required to reproduce the velocity structure of the accreting 
gas within 0.128 pc of the star, providing a better fit than 
the without-disc model. Additionally, as the model with a 
disc was able to adequately reproduce the data, their results 
showed no evidence that ionization or radiation pressure are 
profoundly altering the dynamics of the surrounding accre- 
tion flow. In this work, we carry out the same experiment, 
but instead model the SED and observed infrared images of 
IRAS 20126+4104. 
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In comparison to the molecular line data, the SED and 
infrared images of IRAS 20126+4104 probe different aspects 
of the circumstellar material. The molecular line observa- 
tions are able to measure the velocity of the gas at different 
temperatures, however their density determinations depend 
on abundances that are highly uncertain. In contrast, the 
SED and infrared images probe the dust out to the edge of 
the envelope, which can be used to infer the structure of the 
gas. Hence, comparing the results from both these datasets 
will allow us to obtain a more complete picture of the star 
formation processes occurring in IRAS 20126+4104. 

IRAS 20126+4104 is a well studied example of a nearby 
forming massive star (Lboi = 1.3xlO 4 L , B0. 5 star, |Cesa^| 



roni et al. 1999|). Found within the Cygnus-X star form- 



ing region, it is assumed to be at the same distance as 



the Cygnus OB associations (1.7kpc, |Massey Sz Thompson 
1991; Wilking et al.|1990| ). Several observations of molecular 



tracers at millimetre wavelengths have uncovered a ^0.25 pc 
NW-SE inner jet feeding a large-scale North-South bipolar 



outflow (Cesaroni et al.||1997 


1999 Su et al.||2007| Lebron 


et al. 


2006 


Shepherd et al. 


2000 


). In addition, high res- 



olution centimetre continuum images of the ionized gas to- 
ward IRAS 20126+4104 have been presented by |Hofner et al.| 
( 2007 ). 

Observations of transitions of C 34 S, CH 3 CN, CH 3 OH 



and HCO + (Cesaroni et al. 



1997 



1999 



2005), and NH 3 



( Zhang et al.|1998| KZ10) have uncovered a rotating flat- 



tened disc-like structure, perpendicular to the inner flow, 
with a radius between 800 and 12,000 au, and a gas mass es- 
timated from both dust continuum or molecular line obser- 
vations of 0.65-10 M©. Further evidence of a disc is also pro- 
vided by near and mid-infrared images of IRAS 20126+4104 
(Sridharan et al. 2005 De Buizer 2007), which show a 



dark lane with a similar size and orientation to the disc 
seen in radio observations. A large (radius~1000 au) ro- 
tating structure traced by OH masers was also found by 



Edris et al. (2005), which is coincident with the dark lane 
observed in the infrared. However, recently |de Wit et al.| 
(2009) have presented a high-resolution (0.6") 24.5/xm im- 
age of IRAS 20126+4104, at a slightly longer wavelength 



than those presented in De Buizer (2007 12.5 and 18.3/xm), 
which does not contain a dark lane. 

The infrared to sub-millimetre SED of 
IRAS 20126+4104 has been previously modelled on 



several occasions. Both [Cesaroni et aL ( |1999| ) and Shepherd 



et al. ( 2000 ) modelled the far-IR to millimetre SED using a 
simple greybody to determine the dust temperature and the 
dust emissivity index. Williams et al.| ( 2005| , [de Wit et al 



(2009} and |van der Tak et al.| ( |2Q0Q| ) have also modelled 
the SED of IRAS 20126+4104 using spherically symmetric 
models, to determine the properties of the envelope. More 
recently, iHofner et al.l (120071) modelled IRAS 20126+4104 



using the radiation transfer code previously employed in 
Olmi et al. (2003). Their model consisted of a spherical 



halo with a power-law density between radii of 850 au and 
0.54 pc, plus a constant density cylindrical edge-on disc 
with radii between 34-850 au and a height of 640 au. The 
model provided a good fit to the data down to ~8/xm, 
however, they suggested that adding further components 
to the model, such as bipolar cavities and a more complex 
disc geometry, would improve the fit, especially in the near- 
to mid-IR section of the SED. In this paper, we improve 



upon previous infrared modelling and derive complimentary 
information about the properties of IRAS 20126+4104 not 
previously obtainable from molecular line observations. 

We present the observed SED, infrared images and mea- 
sured brightness profiles of IRAS 20126+4104 in Section [2] 
Section [3] describes the radiation transfer dust code used 
to model them; Section [4] details the modelling, including 
model input assumptions, fitting and the genetic search al- 
gorithm used to search for the best-fitting set of model pa- 
rameters; and Section [5] presents our results, including a 
comparison to the results of KZ10. We give our conclusions 
in Section [6] 



2 THE DATA 

2.1 Near-IR to sub- millimetre SED 

The observed near-IR to sub-millimetre SED of 
IRAS 20126+4104 is shown in Figure [l] The overplot- 
ted solid and dashed lines (respectively the best-fitting 
models for an envelope plus disc and without a disc) are 
discussed further in Section [5] The SED data points were 
collected from the literature; Table [I] lists the wavelength, 
flux and reference for the data displayed in Figure^ If the 
fluxes were derived from apparent magnitudes, these are 
also given in Table [I] We attempted to make sure all fluxes 
contained all of the flux from the source. For instance, 
we used the 2 Micron All Sky Survey (2MASS) extended 
source catalog fluxes for the near-IR. 



2.2 Infrared images 

A three-colour Infrared Array Camera (IRAC) image of the 
mid-IR emission observed toward IRAS 20126+4104, previ- 
ously published by |Qiu et al.| ( [2008] ), is shown in the middle 
panel of Figure [2] where red: 8.0/mi, green: 4.5/mi and blue: 
3.6/xm. In this figure, two lobes of emission oriented simi- 
larly to the NW-SE jet observed bylSu et al.l (120071) can be 



seen, which are most obvious at 3.6/xm and 4.5/xm, while the 
8.0/xm IRAC emission mostly traces the central source. 

To obtain IRAC fluxes for IRAS 20126+4104, we per- 
formed irregular aperture photometry on the images pre- 
viously published in |Qiu et ah] ([2008), using the post- 
Basic Calibrated Data (post-BCD) mosaics produced by the 
Spitzer Science Center (SSC) pipeline version S18.7. The 
photometry was carried out on the short exposure (0.4 s) im- 
ages, as there were several artifacts in the longer exposure 
(10.4 s) mosaics. The apertures were drawn to best avoid 
any foreground stars in the images. 

Measurement uncertainties were <5% at 3.6 and 
4.5 jLtm, ~10% at 5.8 //m, and ^30% at 8 /mi, due to the 
increasing variance of the background in the longer wave- 
length bands. In addition, this should be combined with a 
calibration uncertainty of 10% for pipeline post-BCD mo- 
saic^] Further, due to scattering of flux throughout the im- 
age by the IRAC optics, fluxes of extended sources are highly 
uncertain. Therefore, at a conservative estimate, these fluxes 



1 Section 4.3 of the IRAC Instrument Handbook, version 1.0 
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Table 1. Observed near-IR to sub-millimetre fluxes for IRAS 20126+4104, collected from the literature. 
Apparent magnitudes are given where the flux was derived from these magnitudes. References are given 
in parentheses: (1) |Jarrett et al.| (|2000|), (2 ) thi s work, (3) |Qiu et al.| ( |2008|>, (4) |Price et al.] (2001) , (5) 
|Cesaroni etaT] ( |1999| >, (6) |De Buizer| ( |2007| >, (7) |de Wit et al.| ( |2009^ (8) point Iras Science| ( |1994| > 



Wavelength 


Flux 


Apparent 


Description / Origin 


(/mi) 


(mJy) 


magnitude 




1.235 


4.7 ± 49% 


13.939 ± 0.527 


2MASS Extended Source Catalog, J band (1) 


1.662 


30.1 ± 8.8% 


11.333 db 0.095 


2MASS Extended Source Catalog, H band (1) 


2.159 


134 ± 2.3% 


9.242 ± 0.025 


2MASS Extended Source Catalog, K band (1) 


3.6 


790 




IRAC band 1, photometry: (2), image: (3) 


4.5 


2.1xl0 3 




IRAC band 2, photometry: (2), image: (3) 


5.8 


1.9xl0 3 




IRAC band 3, photometry: (2), image: (3) 


8.0 


1.4xl0 3 




IRAC band 4, photometry: (2), image: (3) 


8.28 


965 ± 4.2% 




MSX band A (4) 


10.2 


320 




UKIRT, MAX (5) 


12.13 


1.10x10 s ± 8.2% 




MSX band C (4) 


12.5 


1.89x10 s ± 10% 




Gemini North, Michelle (6) 


14.65 


Q QQy 1 f|3 _L a nOY 




MSX hand D (A\ 


18.3 


2.35X10 4 ± 15% 




Gemini North, Michelle (6) 


19.9 


3.0xl0 4 




UKIRT, MAX (5) 


21.3 


4.43X10 4 db 6.0% 




MSX band E (4) 


24.5 


6.0xl0 4 




Subaru telescope, COMICS (7) 


25 


1.09x10 s db 5% 




IRAS 25/im (8) 


60 


1.38X10 6 ± 12% 




IRAS 60/im (8) 


100 


1.95x10 s db 13% 




IRAS 100/rni (8) 


350 


2.92xl0 5 ± 9.0xl0 4 




JCMT, SCUBA (5) 


443 


1.62xl0 5 ± 2.4xl0 4 




JCMT, SCUBA (5) 


750 


2.5xl0 4 ± 5xl0 3 




JCMT, SCUBA (5) 


863 


1.9xl0 4 ± 3xl0 3 




JCMT, SCUBA (5) 




1 10 100 1000 

A(/xm) 



Figure 1. The Spectral energy distribution (SED) of IRAS 20126+4104. The best-fitting models for an envelope plus disc and without 
disc (overplotted solid and dashed lines respectively) are discussed in Section [5] The errors shown are those reset to 10% if the error in 
the measured flux was <10%, and the errors in the 2MASS JHK fluxes to a factor of two. The IRAC 4.5/im flux was assumed to be an 
upper limit due to possible contamination by H2 line emission. 
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Figure 2. Left panel: Model IRAC image for envelope plus disc model. Middle panel: Observed IRAC three-colour RGB image of 
IRAS 20126+4104, previously presented by |Qiu et al.| |2008| >. Right panel: Model IRAC image for envelope without disc model. The 
model images have been normalised to the total integrated fluxes given in Table ^ so that the morphology of the emission can be easily 
compared. Stretch: red: 8 /mi, 0-1200 MJysr -1 ; green: 4.5 /im, 0-1800 MJysr -1 , blue: 3.6 /mi, 0-670 MJysr -1 



are accurate to within a factor of two. However, this is ad- 
equate for the purpose of fitting a model to the observed 
SED. 

The middle panels of Figure [3] also show the observed 
K-band, 12.5 /im, 18.3 /im, and 24.5 /im images presented 
in |Sridharan et al.| fl2005| ), |De Buizer| ( |2007| ) and |de Wit| 
et al. (2009), along with corresponding model images (left 
and right panels) which will be discussed in Section [5] 



2.3 Brightness profiles 

Flux profiles of IRAS 20126+4104, summed across a 15.1" 
thick strip aligned with the source outflow axis (PA=303.5 
degrees), were measured for the four IRAC bands and at 
the two wavelengths observed by De Buizer (2007): 12.5 and 
18.3 /im. The normalised profiles are shown in Figure [4] along 
with the best-fitting models to both the SED and profiles, 
for envelope plus disc and without disc models (blue solid 
and red dashed lines respectively), which will be discussed 
further in Section [5] The errors in the profiles shown in each 
panel of Figure [4] reflect the uncertainty in the background 
underlying the source emission due to background fluctu- 
ations, and are calculated as the standard deviation of the 
profiles measured in two strips either side of the main profile, 
which we assumed to contain minimal source flux. 



3 THE RADIATION TRANSFER DUST CODE 

In this paper, we model the SED of IRAS 20126+4104 be- 
tween near-IR and sub-mm wavelengths using the 2D Monte 
Carlo dust radiation transfer code previously employed in 



Wood et al. (2002) and modified in Akeson et al. (2005). For 



a given source and surrounding dust geometry sampled in a 
spherical-polar grid, the code models the anisotropic scat- 
tering and thermal emission by /from dust, calculating radia- 
tive equilibrium dust temperatures (using the technique of 
Bjorkman & Wood 2001), and producing SEDs and multi- 
wavelength images. Dust and gas are assumed to be coupled 



in the code, from which we aim to probe the bulk of the 
material which surrounds IRAS20126+4104. In this section, 
we first describe the density structure of the disc, envelope 
and outflow cavity in the model, then we detail the heating 
of the star via accretion on to its surface, and the dissipation 
of energy via disc accretion. 

We model the circumstellar geometry of 
IRAS 20126+4104 with the same envelope plus disc 
geometry successfully employed to model the SEDs and 



scattered light images of low-mass protostars (e.g. Robitaille 
|et al.|[2007| |Wood et alT][2Q02l jWhitney et~al]|2003| ), and 
also the molecular line emission from IRAS 20126+4104 
(KZ10). The two dimensional flared disc density is described 
between inner and outer radii R m in and iodise by 



Pdisc(^, Z) 



: Po 



exp 



h(zu) 



a) 

where R* is the stellar radius, vj is the cylindrical radius, z 
is the height above the disc midplane and po is set by the to- 
tal disc mass Mdi sc , by integrating the disc density over 2, vj 
and 0, the azimuthal angle. The scaleheight h{vj) increases 
with radius: h{w) — ho {m/R*Y where ho is the scaleheight 
at R* and the disc is assumed to be in hydrostatic equilib- 
rium at the dust sublimation radius. We have assumed the 
parameter /3 = 1.25, derived for irradiated discs around low 
mass stars ( |D'Alessio et al.|1998| , and have taken a — 2.25, 
to provide a surface density of E ~ vj~ x . 

The density of the circumstellar envelope is taken to 
be that of a rotationally flattened collapsing spherical cloud 
flUlrich|[l976l |Terebey et aT|l984| ) with radius Cv X , 

3/2 / it \ -1/2 

penv(r, /i) = . , _^ en ^ 1/0 ( TT- ) ( 1 + 



4vr (GM*Rl) 



1/2 



+ 



P_ 

Mo 



(2) 



where M env is the accretion rate through the envelope, M* 
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Figure 3. Observed (middle panels) and corresponding envelope plus disc and without disc model images (left and right panels respec- 
tively) for several near- and mid-IR observations of IRAS 20126+4104. From top to bottom: K band (Sridharan et al. 2005J, 12.5/im and 
18.3/xm ( pe Buizer|2007| ), and 24.5/mi \de Wit et al.|2009| >. The model images have been normalised to the total integrated fluxes given 
:|1] so that the morphology of the emission can be easily compared. Image stretches (top to bottom): 0-80.7 MJysr -1 , 0-8.98 



in Table _ 
xlO 4 MJysr" 1 



0-2.67X10 5 MJysr" 1 , 0-2.03xl0 6 MJysr~ 



© 0000 RAS, MNRAS 000, 000-000 



6 K. G. Johnston, E. Keto, T. P. Robitaille and K. Wood 



o 

Sh 

a 

x 



CD 



O 




0.8 
0.6 
0.4 
0.2 
0.0 



1 1 1 
12.5 /mi 


1 1 1 

i 


if 






i 

II 


ll 



35 



40 



50 



55 




Offset along profile axis (arcseconds) 



Figure 4. Black error bars: normalised flux profiles for the four IRAC bands and the two wavelengths observed by |De Buizer ( |2007| >: 
12.5 and 18.3 /im. Blue solid and red dashed lines: the profiles of the best-fitting models of both the observed SED and profiles, for 
envelope plus disc and without disc models respectively. 



is the stellar mass, r is the spherical radius, R c is the cen- 
trifugal radius, \i is the cosine of the polar angle (/i = cos#), 
and no is the cosine of the polar angle of a streamline of 
infalling particles in the envelope as r — »• oo. The equation 
for the streamline is given by 



$ + fjio(r/R c - 1) - v(r/R c ) = 



(3) 



which can be solved for /io- In our model, we assume the 
centrifugal radius is also the radius at which the disc forms 



To reproduce the near-IR (Sridharan et al. 2005) and 
IRAC flQiu et al.|2008) images of IRAS 20126+4104, we also 
include a bipolar cavity in the model geometry with density 
Pcav, and a shape described by 

-r^cav 



z(w) — ZqVO 



(4) 



where the shape of the cavity is determined by the parameter 
frcav which we set to be b cav = 1.5, and zo is chosen so that 
the cavity half-opening angle at z = 10, 000 au is C av 

The inner radius of the dust disc and envelope i? m in was 
set to the dust destruction radius R su b , found empirically by 



Whitney et al. (2004) to be 



Rsub — R*(T su h/T*, a 



(5) 



where the temperature at which dust sublimates is assumed 
to be T su b = 1600K, and T*, acc is the effective tempera- 
ture of the star, including heating by accretion shocks on 
to the stellar surface. T*^ cc was found by assuming that 
the material accreting through the disc is in free-fall from 
the inner radius of the gas disc R gas to R* and that half 
of the energy lost by free-fall goes into heating the stellar 



surface (Calvet & Gullbring 1998). The other half of this 



energy is assumed to go into X-ray emission, which is not 
included in the code. As virtually all of the X-rays are re- 
processed for embedded sources like IRAS 20126+4104, the 
determined disc accretion rates may be a factor of two too 
high, as without X-ray emission the disc accretion rate has 
to be a factor of two higher to reproduce the same accretion 
luminosity. We assume .Rgas = 5R* , in keeping with findings 
for discs around low-mass stars ( Shu et al.|1994 ), however we 
note that magnetic disc truncation may not occur for mas- 
sive stars. Therefore, the accretion luminosity which goes 
into heating the stellar surface is 



Lhe 



GM*M dii 



Rcr 



(6) 
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and the effective temperature of the star is then 

T,, acc =T, ( L *^f heat ) 1/4 , (7) 

where L* is the luminosity of the star. 

In addition to luminosity from the star, the dust disc 
emits due to accretion. The energy dissipated through dust 
radiation in a volume element between R m i n and iodise is 



given by the a-disc prescription (e.g. Shakura & Sunyaev 



1973||Kenyon fe Hartmarm||1987| [Whitney et al.||2003| ), 

(8) 

with the wavelength of this emission determined by the 
local dust temperature and opacity. 

For the circumstellar dust properties we adopt the dust 
opacity and scattering properties determined by |Kim et al.| 
( 1994 hereafter KMH), which have a grain size distribution 
with an average particle size slightly larger than the diffuse 
ISM, and a ratio of total-to-selective extinction Ry = 3.6. 

In order to compare the models to existing observations, 
synthetic model images were produced by the dust code 



using the peel- off approach (originally described in Yusef- 
|Zadeh et al.| fl984), whereby each time a photon interacts, 
the probability that the photon reaches the observer, and 
therefore the emergent intensity from that point, is found. 

To produce the model IRAC images, we used the IRAC 
filter spectral response curves, obtained from the SSC web- 
sit^] and convolved the images with the FWHM of the 
IRAC point-spread function (PSF) for each band. 

To produce the images at the four wavelengths shown 
in Figure[3](K band, 12.5/xm, 18.3/xm and 24.5/xm), we ap- 
plied the correct band pass filters, obtained from the United 
Kingdom infrared telescope (UKIRT^ the Gemini telescope 
web page^] and de Wit (private communication), and con- 
volved the images to the observed resolutions (K band: 0.3" , 
12.5/mi: 0.33", 18.3/mi: 0.48" and 24.5/mi: 0.6"). 



4 MODELLING 

In this section, we firstly describe the model input assump- 
tions, including the selection of sensible physical parameter 
ranges for IRAS 20126+4104. Secondly, we detail the \ 2 fit- 
ting procedure carried out to compare the models to the 
data, and finally we outline the genetic search algorithm 
which was used to search parameter space for the most suit- 
able model. 



4.1 Input assumptions: model parameters 

A set of plausible ranges for parameters describing the 
model, in which the genetic algorithm (described in Sec- 
tion 4.3 ) searched for the best fit, is given in Table [5] Ranges 
for the first six parameters - stellar mass, envelope outer ra- 
dius, cavity half-opening angle, inclination, disc mass and 



disc radius - were chosen using results from the literature, 
specifi c ally |Cesaroni et al.| dl997[ [l999| |2005), |De Buizer 



( 2007} , |Edris et al.|fl2005| ), |Qiu et al.|fl2008D and |Zhang et al 
(1998b. Large ranges were assumed for the disc mass and 
radius, as although a disc-like rotating structure has been 
detected at the centre of IRAS 20126+4104 



(e.g. 



Cesaroni 



et al. 2005), it is possible that these observations are also 



tracing part of the envelope. 

The following three parameters in Table[2]- the envelope 
and disc accretion rates and cavity ambient density - were 
also given large but physically plausible ranges. 

The final two parameters, the stellar radius and tem- 
perature, were determined from the sampled stellar mass 
for that model. To do this, we used the M oc M 2 evolution- 
ary track presented in Figure 1 of |Keto fe Wood] j|2006), 



which results from the stellar evolution models of [Chiem 
&; Straniero| ( |1989| ). These calculations show that, before 
taking into account heating from accretion shocks, the pro- 
tostellar temperature and luminosity for stars which reach 
masses greater than 4M© are equivalent to a zero age main- 
sequence (ZAMS) star. This can be seen in Figure 1 of 



Keto 



&; Wood| ( |2006| , where for all the evolutionary tracks (or ac- 
cretion rates), stars which reach masses greater than 4M 
follow the main-sequence line until their accretion is unable 
to counteract the production of helium in their cores, and 
they begin their post-main-sequence evolution. 

Once the stellar mass, radius and temperature were de- 
termined, we calculated the stellar temperature including 
heating from shocks T*, aC c and the surface gravity of the 
star. We then used the nearest solar metallicity [Z — 0.02) 
stellar model atmosphere from the grid of Castelli Ku-| 
|rucz| ( |2004| , or a blackbody spectrum if the temperature 
was larger than 50,000 K, to provide the spectrum of the 
central source. 



4.2 Model fitting 

As part of the genetic algorithm described in Section |4.3[ 
the models were fit to the data using the fitting algorithm 
described in iRobitaille et al. (2007}, where A v , the external 



2 http://ssc.spitzer.caltech.edu/irac/calibrationfiles/spectralresponse 

3 http:/ /irtfweb. ifa.hawaii.edu/~nsfcam/filters. html 

4 http://www.gemini.edu/sciops/instruments/michelle/imaging/filters 



interstellar visual extinction, was a free parameter in the 
fit. We set the distance to be 1.7 kpc, in keeping with the 
adopted distance to IRAS 20126+4104 (Wilking et al.|1990) , 
and allowed the extinction A v to vary between and 40 
magnitudes. 

Before fitting these models to the data, the observed 
flux errors were reset to 10% if they were below this value, 
to account for uncertainties not included in the photomet- 
ric errors such as those due to variability. In addition to 
this, the 2MASS flux errors were reset to a factor of two, 
equivalent to down-weighting these data points by a factor 
of approximately 50 in the % 2 calculation, so that these data 
points would not over-dominate the fit. As shocked H 2 emis- 
sion may be contributing to the measured flux in the 4.5 /am 
IRAC band, we made this value an upper limit by rejecting 
models which had a 4.5 /mi flux higher than the measured 
flux. Each data point was also weighted by the distance in 
log wavelength between its two adjacent data points. This 
was done to allow the algorithm to search for a model which 
reproduced the shape of the entire SED, not only the sec- 
tions with a high density of data. 

The models were also simultaneously fit to the six flux 
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profiles shown in Figure [4] measured from the four IRAC 
images and the 12.5 and 18.3 /mi images in |De Buizer et ah] 
(2005). The model flux profiles were found in the same way 
as the observed profiles, by summing the flux in the model 
images within a 15.1" wide strip aligned with and centred 
on the outflow axis, and normalising them to the total flux. 
When fit, the model profiles were allowed to be shifted by 
±2" (with the same shift for all profiles) around the centre 
of the observed profile to produce the best fit. This was 
necessary since the position of the central source could not 
be accurately determined from the images. 

The combined % 2 for each model given as input to the 
genetic algorithm was calculated as 



2 _ 2 .2 
X — XSED + Xprofiles 



(9) 



where Xsed is the best-fitting reduced \ 2 value for the 
SED alone, and Xprofiles is the overall reduced x 2 for the 
profile fits. 



4.3 The genetic algorithm 

In the context of SED modelling, genetic search algorithms 
(Holland 1975; Goldberg 1989) have previously been used 
by |Hetem &; Gregorio-Hetem| ( [2007| to model the proto- 
planetary discs of T Tauri and Herbig Ae/Be stars. Here we 
describe our algorithm, which was used in conjunction with 
the radiation transfer dust code described in Section [3] to 
search for the best fitting model for IRAS 20126+4104. 

Genetic algorithms are an alternative to traditional op- 
timization methods, such as Markov chain Monte Carlo or 
simulated annealing. They move towards the optimum solu- 
tion in an evolutionary manner, by generating generations 
of solutions or models, based on a given fitness statistic such 
as a x 2 value. Along with simulated annealing, they have the 
advantage that they do not become stuck in local minima, 
as their solution probabilistically tends towards the global 
optimum. In addition, genetic algorithms suit our specific 
problem because they involve producing or running several 
models at the same time. As each of the models was run with 
one million Monte Carlo photons (actually energy packets, 
Bjorkman & Wood 2001), taking approximately 30 minutes 
and 9 minutes per model to complete for the models with 
and without a disc respectively, this allows us to explore pa- 
rameter space much more quickly than optimization codes 
which run one model at a time. 

To create the first generation of models, N models were 
generated by randomly sampling the parameters linearly or 
logarithmically within the ranges given in Table [2] The type 
of sampling (linear or log.) is also noted in the table. Here 
we took N 100. 

The second generation consisted of M models, where 
M was some fraction of N . Here we chose M=20. Models 
in this and subsequent generations were computed by either 
the mutation of one existing model, or the crossover of two. 
We used tournament selection to select which models were 
mutated or crossed. To create a tournament, a fraction of 
the models, k (we chose k = 10%), were selected at random 
from the gene pool of N models. The winner of each tourna- 
ment was determined by first sampling a random number q 
between and 1. The probability p, between 0.5 and 1, which 



we chose to be p = 0.9, was used to determine how likely a 
model is to win the tournament: if q < p, the model with the 
lowest x 2 value was selected. If not, and q < p(p — 1), then 
the second best model was selected, and so on such that if 
q < p(p — l) n_1 , where n is the nth model in the list or- 
dered by x 2 5 the nth model wins the tournament. With this 
probabilistic selection, models that increase \ 2 are occasion- 
ally promoted so that the evolution of the final model is not 
always toward a minimum. Therefore the algorithm allows 
escape from local minima. For mutations, one tournament 
was required; for crossovers, the pair selected for crossing 
resulted from two tournaments. 

A mutation was carried out by selecting one of the var- 
ied parameters at random and resampling it from the origi- 
nal ranges given in Table [2] A crossover was carried out by 
'crossing' the parameters of the two selected models, by ran- 
domly sampling (either linearly or logarithmically) a value 
between the original parameters of the two models. In all 
generations following the first, 75% of the models were the 
result of mutations, and 25% were the result of crossovers. 

When the M models of the next generation were added 
to the gene pool, the M worst (highest % 2 ) models were 
removed so that the gene pool always contained N mod- 
els. All following generations were created using the same 
method as the second. The code was stopped when the best 
fitting model produced a less than 5% reduction in % 2 over 
20 generations. 



5 RESULTS 

The genetic code was run twice - firstly with the model pa- 
rameter ranges given in Table [2] as input for all parameters, 
which we refer to as the 'envelope plus disc' model, and sec- 
ondly with the parameter ranges given in Table [2] for all 
parameters except the disc mass and disc accretion rate, 
Mdisc and Mdi sc , which were set to zero and were therefore 
not treated as model parameters; we refer to this second run 
as the 'envelope without disc' model. The envelope without 
disc model was run to ascertain if the SED and images could 
be adequately produced without a disc, i.e. a simpler model 
which had two fewer parameters. For the envelope without 
disc model, iodise was instead referred to as the centrifugal 
radius of the envelope, R c , as these two parameters are in- 
terchangeable in the models. 

The genetic search algorithm codes were stopped af- 
ter 61 and 48 generations for the envelope plus disc and 
without disc models respectively, when the convergence cri- 
teria were reached. This corresponds to 1300 models or 631 
CPU hours for the envelope plus disc model, and 1040 mod- 
els or 152 CPU hours for the envelope without disc model. 
The resulting SEDs and profiles of these models are shown 
against the data in Figures [I] and [4] and the model images 
are compared to the observed images in Figures [2] and [3] 
The parameters of the best-fitting models are also given in 
Table|3] along with the \ 2 values for the SED-only fit, the fit 
to the image profiles, and the total x 2 for both the SED and 
profiles combined. We have included uncertainties on the 
parameters in Table [3] however these are not formal 1-a er- 
rors given by a A% 2 of 1, as these confidence intervals would 
be unrealistically small for each of the fitted parameters, 
and would in most cases only contain the best fitting model. 
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Table 2. Assumed ranges for model parameters as input to the genetic search algorithm. 



Parameter 


Description 


Value/ Range 


Sampling 


M* 


Stellar mass (Mq) 


5-25 


logarithmic 


- a env 


Envelope outer radius (au) 


15,000 - 120,000 


logarithmic 


q 


Cavity half-opening angle at 10^ au (degrees) 


5-40 


linear 


i 


Inclination w. r. t. the line of sight (degrees) 


75 - 90 


linear 


M disc 


Disc mass (Mq) 


0.1 - 15 


logarithmic 


-^disc 


Disc or centrifugal radius (au) 


100 - 10,000 


logarithmic 


M disc 


Disc accretion rate (Mq yr _1 ) 


10~ 7 - 10~ 2 


logarithmic 


Menv 


Envelope accretion rate (Mq yr _1 ) 


10~ 7 - 10~ 2 


logarithmic 


Pcav 


Cavity ambient density (gem -3 ) 


io- 21 - io- 16 


logarithmic 




Stellar radius 


Determined from stellar mass 




T* 


Stellar temperature 


Determined from stellar mass 





In addition, as the model is only a rough representation of 
reality, this uncertainty would probably underestimate the 
true error. Therefore, we have instead quoted in Table [3] the 
ranges corresponding to models with a reduced x 2 l ess than 
Xbest + 20, and have verified that models in this range still 
provided a reasonable fit to the data. 

5.1 High signal-to-noise runs 

The model SEDs shown in Figure [I] are the result of run- 
ning one million photon models. To obtain higher signal to 
noise in the SED and images, we ran the dust code with 10 
million photons with the final set of parameters produced 
by the genetic code for both envelope plus disc and without 
disc models. The SEDs produced by the 10 million photon 
runs for both models are shown in the left and right panels 
of Figure [5] as solid lines, compared to the SEDs from the 
genetic code one million photon runs, shown as dashed lines. 
The SEDs agree well below ^100/im, but are slightly offset 
at sub- mm wavelengths, due to the temperatures calculated 
by the dust code for the outer cells not being fully converged 
after one million photons, as few Monte Carlo photons have 
interacted in these cells. However, as the difference in SEDs 
is smaller than the uncertainties in the observed fluxes at 
the longest wavelengths, this is unlikely to have a significant 
impact on the results of the genetic algorithm. 

5.2 Parameter \ 2 surfaces 

To understand how well the parameters were determined, we 
plotted for x 2 < 500 the reduced x 2 -values of all the models 
run against each parameter, shown in Figures|6]and[7| As the 
first generation of the genetic code was sampled uniformly 
throughout the chosen ranges, we were therefore able to ad- 
equately sample parameter space and recover an approxima- 
tion of the minimum % 2 surface. Therefore, the solid lines 
in Figures [6] and [7] also show histograms of the minimum x 2 
model in each bin. In addition, the grey shaded areas show 
the ranges of parameter values providing a "good" fit (with 
reduced x 2 ~ Xbest <20) given in Table pi From these Fig- 
ures, we can understand to what extent the observed SED 
and profiles of IRAS 20126+4104 constrain the properties of 
the source. 

For the envelope plus disc model, Figure [5] shows that 
the stellar mass M*, the envelope accretion rate M env , and 



the cavity ambient density p cav are well-determined within 
the chosen ranges, having sharp minimum x 2 surfaces. The 
other parameters are not as tightly constrained, however 
the envelope radius R£^, the disc radius iodise and the disc 
mass Mdisc tend to have better fits at higher values. Also, 
the disc accretion rate Mdisc is not well fit by models with 
values above ~ 10 -4 M0yr -1 . The least well-determined 
parameter is the cavity half-opening angle C av 

For the envelope without disc model, Figure [7] shows 
that the parameters are constrained similarly to those in the 
envelope plus disc model. However, differences include that 
envelope accretion rates below ^lO _4 M yr _1 are strongly 
ruled out, there is a deeper minimum towards higher values 
of the cavity half-opening angle, and no models exist with 
X 2 < 500 for values of the cavity ambient density greater 
than 5xl0 -18 gem -3 . In addition, both the inclination and 
centrifugal radius of the envelope R c are very poorly con- 
strained for the envelope without disc model. 

For the envelope plus disc model, we also include plots 
of the reduced x 2 against each parameter for all the mod- 
els with x 2 < 200 for the SED-only fit (Figure ^ and for 
the profiles-only fit (Figure [9}. This allows us to see more 
clearly whether the SED or profiles are constraining a given 
parameter. For instance, it can be seen that the stellar mass, 
envelope outer radius and disc outer radius are better con- 
strained by the SED, and the cavity opening angle, cavity 
density as well as the viewing angle i are better constrained 
by the profiles. The remaining parameters are similarly con- 
strained by both the SED and profiles. 

5.3 Comparison to the observed SED 

We find that the minimum reduced x 2 f° r the best-fitting 
envelope plus disc model is 31.9, and the reduced x 2 for the 
best-fitting envelope without disc model is 78.4, from which 
we infer that the model with a disc provides a better fit 
to the SED and profiles. For the SED, this can be seen in 
Figure ^ where the envelope without disc model does not 
reproduce the observed near- and mid-IR fluxes. 



5.4 Comparison to the observed profiles and 
images 

Here we compare the model profiles and images to those ob- 
served. Figure [4] presents the model profiles for the envelope 
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Figure 5. Left and right: the SEDs produced by a 10 million photon run and a one million photon run for the envelope plus and without 
disc models respectively. 



Table 3. Parameters of the genetic algorithm best-fitting models 



Parameter 


Description 


Value for envelope 
plus disc model 


Value for envelope 
without disc model 


M* 


Stellar mass (Mq) 


12.7±0.0 


11.8±J;§ 


omax 
^env 


Envelope outer radius (au) 


113,000±iooo o 


79,100±|J™ 


$cav 


Cavity half-opening angle at 10 4 au (degrees) 


29.3d=i2i 


34.2±i;i 


i 


Inclination w. r. t. the line of sight (degrees) 


85.2±!;9 


83.4±6.3 


M disc 


Disc mass (Mq) 


5.90±i;|f 




R disc or R c 


Disc or centrifugal radius (au) 


9200±|0g 


1040±|040 


M disc 


Disc accretion rate (Mq yr _1 ) 


2.58± 4 ;|0xl0- 5 




Menv 


Envelope accretion rate (Mq yr — 1 ) 


3.86±J-gf xlO- 4 


5.35±g-^xl0- 4 


Pcav 


Cavity ambient density (gem -3 ) 


9.55±i 3 4 fxl0- 19 


6.27±g-22 xl0 -i9 


R* 


Stellar radius (R©) 


4.53 


4.34 


T* 


Stellar temperature (K) 


28,700 


27,800 


X 2 (SED) 


X 2 value for the SED fit only 


7.0 


24.6 


X 2 (profiles) 


X 2 value for the fit to the profiles only 


24.9 


53.8 


X 2 (Total) 


Total x 2 for both SED and profiles 


31.9 


78.4 



Note: the uncertainties on the parameter values are given by the range of models with a reduced 

x 2 - xLst < 20 - 



plus disc and without disc models. The left and right pan- 
els of Figure [2] also present three-colour RGB model images 
of the IRAC emission toward IRAS 20126+4104: red: 8 /mi, 
green: 5.8 /im, and blue: 3.6 /xm, for both envelope plus disc 
and without disc models. 

Both the IRAC images in Figure [2] and the profiles in 
Figure [4] show that while both models are able to reproduce 
the observed angular size of the mid-IR emission, the enve- 
lope without disc model is not able to fully reproduce the 
shape of the IRAC emission, as its 8 /im model image is too 
centrally peaked. 

Figure [3] also presents the observed and corresponding 
model images at the wavelengths and spatial resolutions ob- 
served in Srid haran et al.| ( |2QQ5| ) , |De Buizer] ( |20Q7| ) and |de 



Wit et al. 



K band, 12.5/xm, 18.3/im and 24.5/im. All 
model images in Figures [2] and [3] were roughly aligned to the 



observed images using the centre of the dark lane in the ob- 
served and envelope without disc model 12.5/xm images, and 
were rotated clockwise by 56.5° to reproduce the observed 
position angle. 

The morphology of the K band model images (top left 
and right panels of Figure [3| are able to roughly reproduce 
the two outflow cavities imaged by Sridharan et al. (2005), 



shown in the top middle panel of Figure [3] Comparing the 
mid-IR model images shown at the remaining three wave- 
lengths, we see that the mid-IR model images are strongly 
affected by the presence of a disc, specifically in the separa- 
tion of the two outflow cavity emission lobes. In Figure [3] we 
also note that the mid-IR images produced by the envelope 
without disc model reproduce the observed mid-IR images 
much better than those of the model with a disc, and note 
that the lack of a dark lane in the observed 24.5/mi image is 



© 0000 RAS, MNRAS 000, 000-000 



Multi-wavelength modelling of IRAS 20126+4104 11 




Figure 6. Plots of the x 2 -value for models with reduced x 2 < 500 against the nine varied model parameters for the envelope plus disc 
model. The solid line shows a histogram of the minimum % 2 -value in each bin, which provides a rough outline of the minimum x 2 surface. 
The grey shaded areas show the ranges of parameter values providing a "good" fit (with reduced x 2 — X 2 est <20) given in Table [3] 



also reproduced by the envelope without disc model, while a 
dark lane is still seen in the shorter wavelength images. Here 
the shadowing produced by a rotationally flattened envelope 
appears to be adequate to match the observed dark lane at 
12.5 and 18.3/xm. 

However, there is further emission to the north west in 
the observed 12.5, 18.3 and 24.5 /xm images, cut off by an- 
other 'dark lane', which cannot be reproduced by our mod- 
els. Several possibilities which may explain this image mor- 
phology include: 



(i) There is a disc toward IRAS 20126+4104, in agree- 
ment with the observed SED and IRAC images. However, 
the disc is truncated or removed inside the radius delin- 
eated by the dark lane between the central emission and 
the northwest emission (at a radius of ^1.1" or ^1800 au), 
so that the emission within this radius resembles the model 
without a disc. Mechanisms which would truncate the disc 
include stellar winds, ionization and heating, which would 
therefore imply a difference between IRAS 20126+4104 and 
its lower- mass counterparts. 
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Figure 7. Plots of the x 2 -value for models with reduced x 2 < 500 against the seven varied model parameters for the envelope without 
disc model. The solid line shows a histogram of the minimum % 2 -value in each bin, which provides a rough outline of the minimum 
X 2 surface. The grey shaded areas show the ranges of parameter values providing a "good" fit (with reduced x 2 — Xbest <20) given in 
Tabled 



(ii) Precession of the outflow axis with time, so that the 
outflow cavity at smaller radii is not aligned with the out- 
flow cavity at larger radii, explaining the observed disconti- 
nuity in the images at a radius of ^1800 au. In fact, there 
is evidence that the outflow IRAS 20126+4104 is precessing 
( Shepherd et al.|2000| |Cesaroni et"aL|2005 ). 

(iii) The emission within ^1800 au is produced by a dif- 
ferent young star to that producing the larger scale infrared 
emission. A companion is expected, as the binary fraction 
of massive stars is higher than their low-mass counterparts. 



In addition, a binary companion may explain the outflow 
precession. 

5.5 The effect of different dust properties 

To show the effect different dust properties would have on 
the results of the genetic code, and whether different dust 
properties could have a similar effect to the presence of a 
disc, Figure ^] compares the SEDs of the envelope without 
disc model for two different sets of realistic dust properties. 
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Figure 8. Plots of the SED-only -value for models with reduced %L n < 200 against the nine varied model parameters for the 

envelope plus disc model. The solid line shows a histogram of the minimum Xg ED -value in each bin, which provides a rough outline 
of the minimum Xg ED surface. The grey shaded areas show the ranges of parameter values providing a "good" fit (with reduced 

^SED _ ^SED,best <20). 



The solid line represents the SED produced using the same 
dust properties as those used for the genetic code, taken 
from KMH, and the dashed line shows the resultant SED 



using dust properties from Draine (2003a| 2003b[ hereafter 



Draine). The Draine dust model reproduces the Milky- Way 
extinction curve for an i^=5.5 (case A with C/H = 30 ppm) 
using the grain size distribution from Weingartner &; Draine] 
(2001 ). The composition is a mixture of carbonaceous grains 



and amorphous silicate grains, and this model has fewer 



small grains than models that fit a typical interstellar ex- 
tinction curve (such as the KMH model we use in this pa- 
per). This model therefore has a higher albedo at near-IR 



(< 3 fim) wavelengths. Figure 10 shows that the Draine SED 
falls short of the KMH SED in the near- and mid-IR regime, 
and hence a different set of realistic dust properties with a 
higher albedo at near-IR wavelengths is not be able to re- 
produce the effect of the addition of a disc to the geometry 
(such as shown in Figure [I]). This is due to the fact that al- 
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Figure 9. Plots of the profiles-only Xp ro fiies~ vame for models with reduced Xp ro files < ^00 against the nine varied model parameters 
for the envelope plus disc model. The solid line shows a histogram of the minimum Xp ro files~ vame in each bin, which provides a rough 
outline of the minimum Xp ro fiies surface. The grey shaded areas show the ranges of parameter values providing a "good" fit (with reduced 

-^-profiles -^-profiles, best ^20). 



though the increased average grain size increases the albedo 
in the near-IR regime, it also simultaneously increases the 
extinction, which therefore reduces the emergent flux. 

5.6 Comparison with molecular line modelling 
results of KZ10 

In this Section, we compare the properties of 
IRAS 20126+4104 found using two different approaches. 



The first has modelled the profiles of molecular lines ob- 
served toward IRAS 20126+4104 in the millimetre (KZ10), 
and the second has modelled the near-IR to sub-millimetre 
continuum SED and images (this work). As both studies 
find that their models require a disc to adequately fit the 
data, Table [4] compares the envelope plus disc model from 
this work and the parameters of the best-fitting disc model 
of KZ10. 

In Table|4j we see that both models find a similar stellar 
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Figure 10. The SED of the best-fitting envelope without disc model using the dust properties given by KMH (solid line) and Draine 
(dashed line). 



mass, disc mass and disc radius, as well as disc accretion 
rate, which was assumed in the KZIO model to be equal to 
the envelope accretion rate. 

The envelope outer radius was not fit in the KZIO 
model, as their model only extended out to a radius of 
~26,400au or 0.128 pc. In addition, the inclination was as- 
sumed to be 90 degrees, similar to our fitted value of 85.2 
degrees. As in our model, the stellar radius and temperature 
in the KZ10 model were determined from the stellar mass. 

The envelope accretion rate, M env , the envelope den- 
sity at the disc radius, po, and the total mass within a ra- 
dius of 26,400 au, M to t (0.128 pc), all have ^3-5 times larger 
values for our envelope plus disc model than the KZ10 disc 
model. This discrepancy may be explained by the fact that 
our model probes the envelope out to larger radii, and there- 
fore can more accurately determine the density in the enve- 
lope. Alternatively, since the mass density in KZ10 depends 
on an assumed NH3 abundance, this discrepancy may also 
be resolved by increasing the assumed NH 3 abundance by a 
factor of ^3-5. 

We fit the temperature in the midplane calculated by 
the dust code for our best-fitting envelope plus disc model 
with a power-law at radii greater than 10,000 au, and found 
that the temperature varied here as T oc R~ 34 . In com- 
parison, the temperature in the KZ10 disc model was found 
to decrease more steeply than T oc i? _2 ^ 3 , corresponding to 
p < — 1 , where p was used to define the temperature power- 
law exponent in the envelope as T oc R~ 2 ^ 4+p \ 

The disc density ratio, A p is the ratio of the disc den- 
sity to the envelope density at the disc radius (equation 9 
of KZ10). Although the two values of A p given in Table [I] 
differ noticeably, this is due to different values of the enve- 
lope density po, while the disc densities at iodise are similar. 
The disc temperature factor, Bt, is the coefficient in the ex- 
pression describing the disc temperature given in equation 
12 of KZ10. This parameter does not have an equivalent in 
our model, as the temperatures for the envelope plus disc 
geometry were solved for by the dust code. 



Finally, comparing the dynamics of the two models, we 
see that the angular momentum and the velocity at the disc 
radius have similar values for both models. 



5.7 Properties of the disc 

Here we discuss the properties of the disc around 
IRAS 20126+4104 suggested by our results. The disc radius 
^Rdisc for the envelope plus disc model was found by the 
genetic code to be 9200±|§oo au. For a Keplerian velocity 
of 1.11 kms -1 at iodise, the outer radius of the disc would 
therefore take ^3xl0 5 yr to complete one orbit. As this is 
on the order of the time massive stars spend in their accre- 
tion phase, it appears unlikely that the outer regions of the 
disc will have had the time to reach centrifugal or hydro- 
static equilibrium, which occur on approximately the same 
time-scales. However, our results show that, in addition to 
the rotationally flattened envelope, extra mass arranged in 
some form of flattened structure is required, for radii up to 
9200 au, to reproduce the observed SED and profiles. 

As the dust code calculates the temperature of the enve- 
lope and disc geometry, we were able to calculate the Toomre 



Q parameter ( Toomre|l 964) - a measure of stability towards 
fragmentation - in the midplane of the disc as a function of 
radius, 



Q = 



(10) 



where c s is the speed of sound in the gas, Q is the angular 
velocity of the disc, G is the gravitational constant and E is 
the surface density of the disc. We found the surface density 
as a function of radius by integrating equation [I] over the 
disc height z. To find the sound speed, we assumed an ideal 
gas composed of molecular hydrogen, 



kT 



(11) 



2.3m// 

where k is the Boltzmann constant, T is the temperature cal 
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Table 4. Comparison of parameters of genetic best-fitting envelope plus disc model with those from KZ10 



Parameter 


Description 


Value for envelope plus 
disc model (this work) 


Value for model with disc 
(KZ10) 




Stellar mass (Mq) 


12.7+0.0 


10.7 


^env 


Envelope outer radius (au) 


113,000±^oo 


>2o,400 (assumed) 


^cav 


Cavity hall-opening angle at 10* au (degrees) 


29.3+^ 




i 


Inclination w. r. t. the line of sight (degrees) 


or o i 3 Q 


90 (assumed) 


M disc 


Disc mass (Mq) 


5.90±§;§g 


2.5 


-^disc 


Disc or centrifugal radius (au) 


nonnj_300 

9200± 6200 


6900 


M disc 


Disc accretion rate (Mq yr -1 ) 


2.58±5:57 x l° 


(7.6 x 10 -5 ) 


Menv 


Envelope accretion rate (Mq yr -1 ) 


3.86±^-^xl0 -4 


7.6 x 10 -5 


Pcav 


Cavity ambient density (gem °) 


9.55±^4§ 8 xl0 iy 




R* 


Stellar radius (Rq) 


4.53 


4.8 


T* 


Stellar temperature (K) 


28,700 


19,000 


po 


Env. density at Rj (cm -3 ) 


2.39 x 10 5 


7.9 x 10 4 


P 


Temperature power law exponent (T oc i? -2 /( 4 +P)) 


1.8 (>10,000au, T oc R~ - 34 ) 


< -1 (Toe R- 2 / 3 ) 


r 


Angular momentum (au kms ~~ 1 ) 


10,200 


8100 


a p 


Disc density ratio 


0.792 


5.1 


Bt 


Disc temperature factor 




15.0 


Vk 


Velocity at (kms -1 ) 


1.11 


1.2 


Mtot (0.128 pc) 


Total mass within 0.128 pc (M ) 


49.7 


12.6 


Mtot(^ n a v X ) 


Total mass within R^ (M ) 


649 





culated by the dust code and 2.3ra# is the average molecular 
mass of the gas. The Toomre parameter in the disc midplane 
as a function of disc radius for the best-fitting envelope with 
disc model is shown as a black line in Figure [TT) From this, 
we see that the disc is stable (Q > 1), and only reaches a 
value less than unity at radii greater than the disc radius. 
Thus the high temperatures and low densities within the 
disc stabilise it against local fragmentation. We also show 
as grey lines the Toomre parameter for the models with re- 
duced x 2 — Xbest < 20, which, apart from one model, lie 
close to or above 1 at the disc radius. Therefore we conclude 
that across the range of uncertainty in the disc properties, 
the disc around IRAS 20126+4104 is still stable to fragmen- 
tation. 



6 CONCLUSIONS 

In order to determine whether the standard model of low- 
mass star formation can also apply to massive stars, we 
have modelled the near-IR to sub-millimetre SED and sev- 
eral mid-IR images of the accreting embedded massive star 
IRAS 20126+4104, using model geometries which are com- 
monly used to describe forming low-mass stars. 

We have used a Monte Carlo radiative transfer dust 
code to model the continuum absorption, emission and scat- 
tering through two scaled-up azimuthally symmetric dust 
geometries, the first consisting of a rotationally flattened 
envelope with outflow cavities, and the second which also 
includes a flared disc. 

To find the best-fitting set of model parameters, we 
used a genetic algorithm to search parameter space, within 
parameter ranges which bracketed previous results or were 
physically plausible. This also allowed us to produce mini- 
mum x 2 surfaces for each parameter, from which we could 



1000 f — — — — 




0.1 1 — — — — 1 — 1 

1Q _4 10 _ 3 1Q -2 1Q -1 1Q 1Q 1 

Radius in i?disc 

Figure 11. The Toomre parameter in the disc midplane as a 
function of disc radius (in iodise) f° r the best-fitting envelope with 
disc model, shown as a black line, and for the models with reduced 
X 2 - Xbest <20 > snown as g re y lines - 

infer how well the SED and image profiles constrained them. 
We found that the two parameters most constrained by the 
data were the stellar mass and envelope accretion rate, which 
were determined to be ~13 M and ~4x 10 -4 M yr _1 re- 
spectively for the envelope plus disc model, and ^12 M 
and rsj 5 x 10 -4 M yr _1 for the envelope without disc model. 
These two parameters are actually proxies for other more di- 
rectly determined properties of the source - the stellar mass 
in the models is a proxy for the stellar luminosity, and the 
envelope accretion rate is a proxy for the total mass in the 
envelope, which is constrained by the millimetre fluxes. The 
stellar luminosities for both the envelope plus and without 
disc models respectively are 1.3 x 10 4 and 1.0 x 10 4 Lq, 
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and the total masses within R^i^ are 649 and 375 Mq . The 
best-fitting values of the remaining parameters are given in 
Tabled 

Our results show that the envelope plus disc model re- 
produces the observed SED and images better than the en- 
velope without disc model, although the model without a 
disc appears to better-reproduce the morphology of the mid- 
IR emission within a radius of ^1800 au. We have outlined 
several possible causes of this discontinuity, such as inner 
disc truncation or outflow precession. It may be that the 
observed mid-IR images, which cannot be fully explained 
by the envelope plus disc model, are showing the effect of 
the radiation of a massive star on its accreting material. 
However, future observations and modelling are needed to 
determine whether precession and/or a young protostellar 
companion are more likely explanations. 

Comparing our results to those of KZ10, we note that 
both studies find that a model with a disc reproduces the 
observations more successfully than one without. Although 
both studies modelled completely different observations us- 
ing different techniques, we find that our best-fitting model 
parameters are similar to those found by KZ10. While we 
find a higher envelope accretion rate than KZ10, and there- 
fore a higher envelope density and total envelope mass, this 
difference may be due to the fact that the masses determined 
from the molecular line modelling depend on the assumed 
molecular abundances. For example, increasing the assumed 
molecular abundances by a factor of 3-5 would bring these 
values into agreement. In addition, as the observed SED is 
more sensitive to emission from the entire envelope, not only 
the material within a radius of 26,400 au modelled by KZ10, 
we were able to probe the mass in the envelope more accu- 
rately. 

Our best fitting envelope plus disc model has a disc ra- 
dius of 9200 au. As at this radius the disc requires ^3x 10 5 yr 
to complete one orbit, we find it is unlikely that the outer 
regions of the disc have had time to reach centrifugal or 
hydrostatic equilibrium. However, by using the temperature 
along the disc midplane found by the dust code, we calculate 
that this disc would be stable to local fragmentation. 
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